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METHOD OF DETECTING SYSTEM 
FUNCTION BY MEASURING FREQUENCY 
RESPONSE 

CROSS-REFERENCE TO RELATED 5 

APPLICATIONS 

This application is a continuation-in-part of U.S. patent 
application Ser. No. 1 1/825,629, filed Jul. 5, 2007, now U.S. 
Pat. No. 7,395,163, issued Jul. 1, 2008, which is a continua- 10 
tion of U.S. patent application Ser. No. 1 1/3 13,546, filed Dec. 

20, 2005, now abandoned, which claims the benefit of U.S. 
Provisional Patent Application Nos. 60/637,969, filed Dec. 

20, 2004, and 60/724,631, filed Oct. 7, 2005. The disclosure 
of each of these applications are hereby incorporated by ref- 15 
erence in their entirety, including all figures, tables and draw- 
ings. 

STATEMENT REGARDING FEDERALLY 
SPONSORED RESEARCH OR DEVELOPMENT 20 

The subject invention was made with government support 
under a research project supported by NASA, Grant No. 
NNA05AC24C. The government has certain rights in this 
invention. Further, the United States Government has certain 25 
rights in this invention pursuant to Contract No. DE-AC07- 
05ID1 45 1 7 between the United States Department of Energy 
and Battelle Energy Alliance, LLC. 

BACKGROUND OF THE INVENTION 30 

Electrochemical Impedance Measurement Systems use the 
Bode analysis technique to characterize an impedance of an 
electrochemical process. It is a well-established and proven 
technique. A battery being evaluated is excited with a current 35 
that is a single frequency and its response is measured. The 
process is repeated over a range of frequencies of interest 
until the spectrum of the impedance is obtained. The method 
is effective but time consuming, as the process is serial. A 
parallel approach using band width limited noise as an exci- 40 
tation current can obtain the same information in less time. 
The system response to the noise is processed via correlation 
and Fast Fourier Transform (FFT) algorithms and many such 
responses are averaged. The result is the spectrum of response 
over the desired frequency range. The averaging of many 45 
responses also makes this process somewhat serial. Another 
technique assembles the current noise waveform from a sum 
of sinusoids, each at a different frequency. The system 
response as a time record is acquired and processed with the 
FFT algorithm. To reduce noise, multiple time records of 50 
waveforms are processed and their resultant spectra averaged. 
This process is also serial. 

There remains a need for real-time acquisition of battery 
impedance for control and diagnostics over a limited fre- 
quency range. This method of acquisition should be a true 55 
parallel approach that uses a single time record of battery 
response with a duration compatible with a real-time control 
process. 

BRIEF SUMMARY OF THE INVENTION 60 

The invention involves using a parallel approach to analyze 
battery impedance or other system functions. A number of 
frequencies are selected over which the battery is to be tested. 
These frequencies are assembled into an Excitation Time 65 
Record (ETR) that is the Sum of the Sinusoids (SOS) of the 
frequencies and the length of such periods of the lowest of the 


2 

frequencies. The ETR is conditioned to be compatible with 
the battery. The battery is then excited with the ETR and a 
Response Time Record (RTR) is captured. The RTR is then 
synchronized to the ETR and processed by a series of equa- 
tions to obtain frequency response. 

In one preferred embodiment, the RTR is processed to 
obtain estimated frequency components of magnitude and 
phase for one of the selected frequencies. Processing is 
repeated to obtain estimated frequency components for each 
selected frequency. Frequency components are reassembled 
to obtain an Estimated Time Record (ETR). The ETR is 
subtracted from the captured RTR to get an error. The error is 
minimized to achieve the frequency response estimate. Error 
is minimized using Compensated Synchronous Detection 
(CSD) using a CSD algorithm, which can be implemented by 
a neural network. 

In another preferred embodiment, all excitation frequen- 
cies of the SOS are harmonics by powers of two. The sample 
period likewise is a power of two with all the SOS frequen- 
cies. The RTR is rectified relative to a square wave and a 
90-degree shifted square wave of one of the SOS frequencies. 
Integrating the processed RTR results in an “in phase” and 
“quadrature” sum that is easily processed to yield the magni- 
tude and phase shift of the desired frequency components. 
Frequency components are assembled to obtain frequency 
response. 

The subject method allows a parallel implementation for 
swept frequency measurements to be made utilizing a com- 
posite signal of a single time record that greatly reduces 
testing time without a significant loss of accuracy. 

BRIEF DESCRIPTION OF THE SEVERAL 
VIEWS OF THE DRAWINGS 

FIG. 1 shows a Filter Ideal Magnitude Response. 

FIG. 2 shows a Filter Ideal Phase Response. 

FIG. 3 shows a Filter Uncompensated Synchronous 
Detected Magnitude Response. 

FIG. 4 shows a Filter Compensated Synchronous Detec- 
tion (CSD) Magnitude Response. 

FIG. 5 shows a Filter CSD Phase Response. 

FIG. 6 shows a Lumped Parameter Model (LPM). 

FIG. 7 shows a portion of a Sum of Sines (SOS) signal to 
LPM, 13 lines, 10 periods. 

FIG. 8 shows a portion of an LPM time response, 13 lines, 
10 periods. 

FIG. 9 shows an LPM ideal magnitude response, 13 lines, 
10 periods. 

FIG. 10 shows an LPM ideal phase response, 13 lines, 10 
periods. 

FIG. 11 shows an LPM uncompensated magnitude 
response, 13 lines, 10 periods. 

FIG. 12 shows an LPM CSD magnitude response, 1 3 lines, 
10 periods. 

FIG. 13 shows an LPM CSD phase response, 13 lines, 10 
periods. 

FIG. 14 shows a portion of an LPM SOS signal, 25 lines, 10 
periods. 

FIG. 15 shows a portion of an LPM time response, 25 lines, 
10 periods. 

FIG. 16 shows an LPM ideal magnitude response, 25 lines, 
10 periods. 

FIG. 17 shows an LPM ideal phase response, 25 lines, 10 
periods. 

FIG. 18 shows an LPM uncompensated magnitude 
response, 25 lines, 10 periods. 



US 8,150,643 B1 


3 

FIG. 19 shows an LPM CSD magnitude response, 25 lines, 
10 periods. 

FIG. 20 shows an LPM CSD phase response, 25 lines, 10 
periods. 

FIG. 21 shows a Mean Squared Error (MSE) comparison 
for a low-pass filter frequency response. 

FIG. 22 shows an MSE comparison for detection of LPM 
impedance. 

FIG. 23 is a flow chart showing a preferred method of the 
invention. 

FIG. 24A shows a Fast Summation Transformation (FST) 
impedance spectrum as magnitude vs. frequency. 

FIG. 24B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 24C shows an FST impedance spectrum as a Nyquist 
plot of imaginary component vs. real component. 

FIG. 25A shows an FST impedance spectrum as magnitude 
vs. frequency. 

FIG. 25B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 25C shows an FST impedance spectrum as a Nyquist 
plot of imaginary component vs. real component. 

FIG. 26A shows an FST impedance spectrum as magnitude 
vs. frequency. 

FIG. 26B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 26C shows an FST impedance spectrum as a Nyquist 
plot of imaginary component vs. real component. 

FIG. 27A shows an FST impedance spectrum as magnitude 
vs. frequency. 

FIG. 27B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 27C shows an FST impedance spectrum as the 
Nyquist plot of imaginary component vs. real component. 

FIG. 28A shows an FST impedance spectrum as magnitude 
vs. frequency. 

FIG. 28B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 28C shows an FST impedance spectrum as a Nyquist 
plot of imaginary component vs. real component. 

FIG. 29A shows an FST impedance spectrum as magnitude 
vs. frequency. 

FIG. 29B shows an FST impedance spectrum as phase vs. 
frequency. 

FIG. 29C shows an FST impedance spectrum as a Nyquist 
plot of imaginary component vs. real component. 

FIG. 30 is a flow chart showing another preferred method 
of the subject invention. 

DETAILED DESCRIPTION OF THE INVENTION 

The method of the subject invention allows for real time 
estimation of a battery’s impedance spectrum. The shift of a 
battery’s impedance spectrum strongly correlates to the 
health of the battery. Therefore, the subject method provides 
in- situ diagnostics for state-of-health estimation of the bat- 
tery, which is critical for enhancing the overall application’s 
reliability. The subject method measures a frequency 
response of a unit under test, for example, a battery. The 
battery under test is excited by the sum of sinusoids of a 
number of test frequencies. A response time record is cap- 
tured, then processed, to obtain estimates of frequency com- 
ponents for each of the number of frequencies. Estimated 
frequency components are assembled to achieve a frequency 
response. 

A system function can be identified over a limited number 
of specific frequencies at step 10 as shown in the flow chart of 


4 

FIG. 23. The desired frequencies are assembled as an excita- 
tion time record 12 that is a sum of those sinusoids and have 
a length of several periods of the lowest frequency. The time 
step selected must be compatible with Shannon’s sampling 
5 constraints for the highest frequency component. The indi- 
vidual waveforms should be sine waves of equal amplitude 
but alternating signs for a phase shift of 1 80 degrees between 
the components. Alternating the 180-degree phase shift will 
minimize a start-up transient. The Root Mean Square (RMS) 
to and the rogue wave peak (sum of the absolute values of all 
component peaks) of the assembled time record must be 
compatible with the system being excited and the Data Acqui- 
sition System (DAS) that will capture the response. 

The Excitation Time Record (ETR), as a current for imped- 
15 ance identification or voltage for system function identifica- 
tion, is signal conditioned 14 to be compatible with the Unit 
Under Test (UUT). As part of the signal conditioning, anti- 
aliasing filters ensure that all frequencies generated by the 
digital -to -analog conversion process other than the intended 
20 frequencies are suppressed. The UUT is excited 16 by the 
ETR and a time record of the UUT response is captured by the 
DAS. The UUT Response Time Record (URTR) is synchro- 
nized to, and the same length as, the ETR. 

A preferred embodiment of the processing of the response 
25 time record is described below. In order to fit steady-state 
sinusoidal response assumptions, a preselected number of 
data points R at the beginning of the URTR must be discarded. 
In general, the sum of those data points total to a time that is 
larger than the transient response time of the UUT at the front 
30 end of the ETR. The UUT corrected time response is referred 
to as the CTR. 

In order to fit steady-state sinusoidal response assump- 
tions, a preselected number of data points, R, at the beginning 
of the URTR must be discarded. In general, the sum of those 
35 data points total to a time that is larger than the transient 
response time of the UUT at the front end of the ETR. The 
UUT Corrected Time Response is referred to as the CTR. 

The first estimate of components, magnitude, and phase of 
the frequency response is made by processing 18 the URTR 
40 via Equations 3 through 1 0 . It is important that a zero mean be 
established prior to processing with Equations 3 to 10 (see 
below). 

The core of this whole concept is that an estimate of the 
UUT Corrected Time Response, the CTR, is made by reas- 
45 sembling 20 the CTR using the estimates of the individual 
frequency components with the same time step and then 
discarding the first R time steps to become the ECTR. The 
difference between the CTR and the ECTR is an error 22 and 
minimizing 24 this error will increase the accuracy of the 
50 frequency response estimates. 

A first approach to minimizing the error between the CTR 
and the ECTR is Compensated Synchronous Detection 
(CSD). The CSD algorithm synthesizes a residual time record 
of the original time record using the magnitudes of the in- 
55 phase and quadrature components for each frequency, except 
the one to be detected. This synthesized residual is then sub- 
tracted from the original time record. The resulting compen- 
sated time record is processed with synchronous detection 
and a new compensated estimate of the response at the detec - 
60 tion frequency is obtained. Since all of the other components 
in this compensated time record are suppressed, the error 
from leakage at those other frequencies is less. This process is 
repeated for each of the frequencies. Assembling the residual 
time record and generating the compensated time record are 
65 illustrated by Equations 1 1 and 12 (see below). 

Another approach to minimize the error between the CTR 
and the ECTR is to use a neural network. A first estimate of 
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the component magnitudes and phases is made as described 
for the CSD technique. Those values are stored and the ECTR 
is calculated. This signal is then subtracted from the CTR to 
produce a response residual. The synchronous detection is 
then performed upon this residual and the component mag- 5 
nitudes and phases are again stored. These components are 
then used to reconstruct an estimate of the residual signal. 
This estimated residual signal is subtracted from the initial 
residual signal to produce a residual, residual signal. This is 
then synchronously detected and the loop starts again. This is 
repeated as many times as desired, each time with the result- 
ant components being stored. The assumption is that there is 
a functional relationship between these resultant components 
and a true system response. A neural network is then used to 15 
determine this relationship. The previously stored results 
become the test dataset for the neural network. The network 
has been trained previously on a similar unit and a known 
ideal response (e.g., battery impedance measured using elec- 
trochemical impedance spectroscopy). The output of the net- 20 
work is the estimate of the response. 

A further step of the subject method is to shift the complete 
set of desired frequencies and to repeat the whole process. 
This step could be repeated many times with different shifts to 
develop a high-resolution frequency response that for battery 25 
impedance could be comparable to that provided by electro- 
chemical impedance spectroscopy. Thus, the subject method 
can provide capability for both limited frequency response in 
real-time or high-resolution frequency response not in real- 
time for periodic system in-depth diagnostics. 30 

The system of the subject invention is based on the follow- 
ing theoretical design. The Unit Under Test (UUT) is excited 
with a limited sum of sinusoids, each at a different frequency 
that is spread over the range of interest. The magnitude, 
frequency and phase of each sinusoid making up the sum are 35 
known. If a total response of the system is measured via a 
sample data system at an acceptable sample rate and an 
adequate duration time record is acquired, then a simple 
algorithm that uses the known magnitude, frequency and 
phase of each individual sinusoid will process the single time 40 
record. This analysis will obtain the true Bode response at the 
selected frequencies spread over the range of interest all in 
parallel. The following synchronous detection analysis is the 
basis of this simple algorithm. The reference waveform is 
chosen as a sine, as at time zero everything will be at zero. 45 
Equation 1 gives the relationship for a parallel excitation. 


f n (r) = ^ /\ / sin(dt>/ r + (pin,) 


( 1 ) 


50 


Equation 2 gives the measured sampled data response of 
the system 


55 


fom [J] = J] B, sin(w ; ( j - 1 )Ar + (pout -, ) ; j=l:N 


(2) 


60 

Where: 

A z is the amplitude of the \ th input sinusoid 

B. is the amplitude response of the \ th output sinusoid 

oo z is the radian frequency of the \ th sinusoid 

At is the time step of data system 65 

(|)in z is the phase of the \ th input sinusoid 

(|)out z is the phase response of the \ th output sinusoid 


N is the number of points of the response time record 
M is the number of different sinusoids of the excitation 
time record 

Each component magnitude and phase of the system 
response at all the excitation frequencies can be obtained via 
the following synchronous detection analysis. This analysis 
quantifies the response at the k th radian frequency with the 

“in phase” and “quadrature” response. The analysis incorpo- 
rates a feature of discarding a preselected number of points R 
at the beginning of the system response in order to meet the 
assumption of steady-state sinusoidal response. Additionally, 
for most applications, prior to processing the data, the mean 
of the acquired time record should be computed and deleted. 
The presence of a non-zero mean could corrupt an estimate of 
the lowest frequency component. 

In phase: 



A K sin((x) K (j - l)Ar + 
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A k Bk (5) 

F out K = — y— cos {(pin K - (pout K ) 

If the time record were infinite, the summation over j would 
average everything to zero except the final result of Equation 
5. The quadrature analysis follows in the same way. 
Quadrature: 




™ ( 
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out K N-R 
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i*k = 1 

sin((w* - (x)i){j - l)Ar + (pin K - 


Fq o t 


ArBr 

: 2 sin ((pout K - (pin K ) 


( 8 ) 


Again, the summation over j for infinite time record aver- 
ages everything to zero except the final result of Equation 8 . 
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Equations 5 and 8 can be combined to give magnitude and 
phase for the k th frequency response. 


I F out ^ I — -J foutfc + fq ou t 


(9) 


A K B K I I ~ ArBr 

y sin L {(pin K - (pout K ) + cos 2 ((pin K - (pout K ) = — 


between a completely synthesized time record response and 
the original time record is minimized. 

The following examples are offered to illustrate the method 
of the subject invention and should not be construed as lim- 
iting. 


Example 1 


-.(^L tan ' 1 
l tout K ) 


AkBk 

2 sin {(poutK ~ <pin K ) 


AkBk 


cos {(poutK ~ (pinfc) 


( 10 ) 


sin {(f>out K -(f)in K )\ 

11 7T— : tt — - = ($out K - (pin K ) 

l cos ((pout K - (pin K ) ) 


10 


15 


Equations 9 and 10 give the final results for synchronous 
detection. This process is repeated for all M of the excitation 
sinusoids. As stated, the results depend on the time record 
being infinite in duration. This is not the case and, thus, the 20 
results are in error by leakage from the other components. 
This error can be reduced without increasing the time record 
using the following algorithm for Compensated Synchronous 
Detection (CSD). 

The CSD algorithm synthesizes a residual time record of 25 
the original time record using the magnitudes of the in-phase 
and quadrature components for each frequency except the one 
to be detected. This synthesized residual is then subtracted 
from the original time record. The resulting compensated 3Q 
time record is then processed with synchronous detection, 
and a new compensated estimate of the response at the detec- 
tion frequency is obtained. Since all of the other components 
in this compensated time record are suppressed, the error 
from leakage at those other frequencies will be less. This 35 
process is repeated for each of the frequencies. Assembling a 
residual time record and generating the compensated time 
record are illustrated by Equations 1 1 and 12. 


M 

Irk [j] = ( F P sin (u p (j - l)Ar) + F Qp cos((o p (j - l)Ar)); 

p=l,p±K 


40 

( 11 ) 


Analytical Testing on a Sum of Sines 

The CSD algorithm was evaluated using a simple signal 
that was assembled from a finite sum of equal amplitude sine 
waves (Sum of Sines (SOS)) with frequencies distributed 
logarithmically over a limited range. The objective of the 
analysis was to assess how well the CSD algorithm could pick 
out an amplitude for each component. 

To check out the concept analytically, a MATLAB® matrix 
calculation computer software code was written that was a 
logarithmic mix of 5 equal unity amplitude frequencies (5 0 ' 5 , 
5 1 , 5 1 ' 5 , 5 2 , 5 2 5 Hz). The acquired time record was set to 10 
periods of the lowest frequency and the time step was set to 
V 10 of the period of the highest frequency. As per Equation 9, 
error-free detection should estimate the amplitude of 0.5 for 
each component. Table 1 gives an estimate for a first pass or 
simple synchronous detection and a second pass, the CSD 
algorithm. The MATLAB® matrix calculation computer 
software code for this analysis is disclosed by Morrison et al., 
2005, Algorithms as MATLAB Code for Real time Estima- 
tion of Battery Impedance. 


TABLE 1 


Compensation algorithm analytical results 

Frequency 


5.5 5 1 51.5 5 2 

5 25 

1st Pass (Simple Synchronous Detection) 

Amplitude 

0.5060 0.5060 0.4975 0.4988 

2nd Pass (CSD Algorithm) 

0.5008 

Amplitude 

0.5004 0.4998 0.5004 0.5000 

0.5003 


j = R + l:N 


CfKoutDMoutDJ-fRKD] (12) 

Where: 

f out is the original time record 
^RK is the correction time record 
Ci K out i s ^e compensated time record 

is the estimated in phase amplitude response at the p th 
frequency 

F Qp is the estimated quadrature amplitude response at the 
p* frequency 

C 0 p is the radian frequency of the p th sinusoid 
At is the time step of the data system 
N is the number of points of the output time record 
M is the number of different sinusoids of the excitation 
function 

R is the number of points of the output time record that are 
discarded 

The synchronous detection algorithm described by Equa- 
tions 1 through 8 is applied to the compensated time record of 
Equation 12 and better estimates of F^ and F^ are obtained. 
This process can be repeated again until the total difference 


As seen in Table 1 analytically, the compensation tech- 
nique does appear to work as the error for the second pass is 
much reduced. This initial check of the concept was applied to 
detect amplitude only and not phase detection. The signal is a 
sum of equal amplitude sine waves being decomposed into 
the individual components by the algorithm. 

Example 2 

Analytical Testing of a Low-Pass Filter 

A recursive model of a second order low -pass function was 
excited with an SOS input signal. The CSD algorithm was 
then used to estimate the frequency response at each of the 
specific frequencies making up the SOS. 

A spread of 13 specific frequencies was chosen that were 
spaced in % decade steps starting from 0.1 Hz up to 100 Hz. 
Using these frequencies, a mix of equal unity amplitude sine 
waves was created. This range of frequencies was picked as 
research performed at the Idaho National Laboratory with 
batteries is typically over this frequency spread. The signal 
was discretized with a time step that was 10 % of the period of 
the highest frequency. The length of the time record was set at 
10 periods of the lowest frequency. 


45 


50 


55 


60 


65 
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A recursive model of a second order Butterworth low-pass 
function was developed. A center frequency was set at the 
middle of the SOS frequency spread. A filter response to the 
SOS input time profile was computed. 

Finally, the CSD code was developed to estimate the filter 5 
frequency response from the time profile generated by the 
recursive model code. That code is the implementation of 
Equations 1 through 12. Additionally, the implementation has 
the ability to discard a number of user-selected points at the 
beginning of the time profile such that the remaining data to 
better fits the assumption of “steady-state” response. The 
following 6 figures are MATLAB® matrix calculation com- 
puter software discrete plots of an ideal frequency response, 
an uncompensated frequency response, and finally, a com- 
pensated frequency response. All 3 have a magnitude plotted 15 
against a Log of frequency. FIGS. 1 and 2 are an ideal mag- 
nitude and phase response. FIG. 3 is an uncompensated mag- 
nitude response, while FIG. 4 is a compensated magnitude 
response. The improvement seen by comparing the last two 
figures is clear FIG. 5 is the compensated phase response. 20 
Error at the higher frequencies for the compensated result is 
likely due to processing a small signal at the high frequencies 
relative to larger signals at the lower frequencies. The net 
result is that the one-time record yields a limited number of 
points for both magnitude and phase. This technique shows 25 
promise for real-time applications. The MATLAB® matrix 
calculation computer software code for this analysis is dis- 
closed by Morrison et al., 2005, Algorithms as MATLAB 
Code for Real time Estimation of Battery Impedance. The 
next approach will apply the concept analytically in an 30 
attempt to identify the impedance of a computer model of a 
battery. 


Example 3 

Analytical Testing of CSD with a Battery Model 


35 


The CSD algorithm was evaluated analytically via a com- 
puter simulation of the detection of the impedance of the 
Lumped Parameter Model of a battery (LPM) that was devel- 40 
oped by the Idaho National Laboratory (INL) (see, Freedom- 
CAR Battery Test Manual, 2003). A computer model for the 
LPM that will simulate battery voltage response to an arbi- 
trary battery current was also developed at INL by Fenton et 
al . , 2005 . The voltage response of the model normalized to the 45 
current in the frequency domain will be the battery imped- 
ance. The equivalent circuit for the LPM with parameter 
identification is shown in FIG. 6. 

The LPM was excited with a current source I /Ar that was an 
SOS, and the CSD algorithm was used to identify the imped- 50 
ance seen looking into the LPM over a limited range of 
discrete frequencies. It should be noted that the polarity of the 
voltage response was defined as negative because the SOS 
excitation current was a discharge (negative relative to imped- 
ance). The CSD algorithm was used to obtain the frequency 55 
response of the LPM. The CSD response magnitude and 
phase were compared to the ideal response. Initially, the 
algorithm failed to match ideal impedance. The response of a 
battery terminal voltage to a discharge SOS current signal 
will contain a DC term caused by the DC battery voltage. 60 
Synchronous detection of any specific frequency in the 
response will cause a noise frequency in the resultant spec- 
trum at the frequency being detected, at an amplitude of the 
product of the DC battery voltage, and an amplitude of the 
detection signal. That noise amplitude will be large relative to 65 
the signal being detected. Averaging enough cycles in the 
resultant time record will reject this noise. However, for a 


real-time application, the length of the time record needs to be 
short and not long. The problem was resolved when the mean 
was deleted from the prediction response of the LPM. The 
number of frequency lines was set to 13 and were logarith- 
mically spread from 0.01 Hz to 1 .0 Hz. The time record was 
set to 10 periods of the lowest frequency. In the CSD algo- 
rithm, no data points were discarded. Table 2 gives analysis 
specifics with LPM data that is typical for some lithium 
batteries that INL had recently tested. INL performed the 
testing per methods in the FreedomCAR Battery Test 
Manual, 2003. 


TABLE 2 


Representative LPM and Analysis Data 

Voc = 3.8 
Cp= 666.6667 

At Rest 

Coc= 1.6667e+003 

At Rest 

Ro = 0.0250 
Rp= .0150 
M= 13 

Number of frequency lines 

Dt = .01 

Time step, sec 

N = 100000 

Total number of points 

F = .01 

Starting Frequency, Hz 

FF = 10 

Frequency spread in decades 

S = .125 

Step size (log) over the decades, 

NN= 10 

8 steps per decade 

Length of time record in number 


of periods of lowest frequency 


The following 7 figures are the plots of the analysis results. 
FIG. 7 is the time record of the SOS current signal. FIG. 8 is 
the time record of the LPM voltage response. FIGS. 9 and 10 
are the LPM ideal impedance magnitude and phase. FIGS. 11 
and 12 are the uncompensated and the compensated magni- 
tude response. FIG. 13 is the compensated phase response. 

The number of frequency lines was increased to 25 and 
everything else was left the same. The following 7 figures 
illustrate the results. FIG. 14 is the time record of the SOS 
current signal. FIG. 15 is the time record of the LPM voltage 
response. FIGS. 16 and 17 are the LPM ideal impedance 
magnitude and phase. FIGS. 18 and 19 are the uncompen- 
sated and the compensated magnitude response. FIG. 20 is 
the compensated phase response. 

Additional cases run showed that as the length of the 
acquired time record in the number of periods of the lowest 
frequency gets cut back, the number of frequency lines that 
can be resolved without a big increase in error must be cut 
back. For example, 5 periods with 25 lines gave terrible 
results, but 5 periods with 13 lines was fine. These positive 
results are only analytical. Nevertheless, they offer promise of 
positive expected performance when applied to a physical 
system. All the MATLAB® matrix calculation computer 
codes for this analysis are given by Morrison, 2005, Algo- 
rithms as MATLAB® Code for Real time Estimation of Bat- 
tery Impedance. 


Example 4 

Neural Network Enhanced Synchronous Detection 

In order to improve the accuracy of the CSD, studies were 
conducted upon neural network enhancement of the detection 
of the individual frequency components of the response of a 
linear system to an SOS input signal. The concept is very 
similar to the CSD, with some slight changes and a neural 
network output layer. For a second order low -pass filter, ordi- 
nary synchronous detection performed on the filter response 
showed a mean squared error (MSE) of 2.6xl0 -3 . The com- 
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pensated synchronous detection technique displayed an MSE 
of 1.6xl0 -3 . The neural network enhanced synchronous 
detection showed an MSE of 0.2x1 0 -3 . Results for the 
lumped parameter model of a lithium ion battery were similar. 

The theory of neural network enhanced synchronous 5 
detection is based on a classical synchronous detection and an 
inherent error associated with time records of finite length. 
Given an input signal comprised of a sum of N sinusoids: 


12 

estimate would be exact if the time record was infinite, but 
since it is of finite length, the reconstructed signal is not 
exactly correct. This may be viewed as containing noise. 


N 

y[i] = ^ ^„sin(w„/Ar) + y^cos^/Ar) + rj[i] 

n= 1 


10 

N 

x{t) = ^ x{OJ n t) 

n = 1 


the output of the system would be: 


15 


N 

y(0 = 

n = 1 


20 


Being that the input and output are sinusoids, they are 
assumed to have started at t=-oo and continues to t=oo. In 
reality, however, this is not the case. First, the time record of 
the signal is finite in length; second, it is sampled. Given a 25 
sampling frequency of at least the Nyquist frequency (twice 
the highest frequency in the signal), the signal can be recon- 
structed without error. This does not, however, rectify the 
finite time length of the signal. Because of this, errors enter 
into the synchronously detected frequency components. 30 
Synchronous detection involves multiplying the acquired 
signal by sines and cosines of the desired frequencies and 
summing the results, as shown below: 


a u> n = linw y[j]sm[aj n jAt] 

j=~i 

(\o n = lim^co j.j] y[j]cos[oj n jAt] 


35 


40 


The magnitude of a given frequency is obtained as follows: 
M^=?Ja^ 2 +^ 2 45 

and the phase may also be obtained as follows: 



50 


where: 

y[j]=the sampled signal 

a=the sine component of the response for the frequency co„ 55 
p=the cosine component of the response for the frequency co„ 

At=the sample time step 

Notice that the summations are infinite in length. In appli- 
cation, the summation would be from 0 to the length of the 
recorded signal . If the time record was infinite, then any errors 60 
would cancel out. Since our time record is not infinite, errors 
remain in the calculated response. 

It is important to note that an estimate of the original signal 
can be reconstructed by summing the sine and cosine signals 
of the different frequencies, where each sine signal is multi- 65 
plied by its sine response component and each cosine signal is 
multiplied by its cosine response component. The resulting 


where: 

r|[i]=the noise due to error 

At=the sample time step 

i=the \ th position in the acquired time record 

N=the number of frequency lines 

In the compensated synchronous detection approach, this 
difference between the actual signal and the reconstructed 
signal is exploited to filter out all but one of the frequencies 
and to increase the accuracy of the measurement (see above). 
In Neural Network Enhanced Synchronous Detection 
(NNESD), the approach is slightly different. The premise is 
that the residual left from subtracting the reconstructed signal 
from the actual signal still contains some useful information. 

r[i]=y[i]-m 

For each frequency of interest, the sine and cosine compo- 
nents are detected on the original signal. These are then used 
to reconstruct the signal and then stored. The reconstructed 
signal is then subtracted from the original signal to obtain the 
residual signal. Synchronous detection is then performed 
upon this residual signal; once again the sine and cosine 
components are used to reconstruct the residual signal and are 
also stored. This second residual signal is then subtracted 
from the first residual signal. This loop continues until a 
sufficient number (M) of sine and cosine responses are 
obtained for each frequency. Now, the assumption that is 
made is that there is some functional relationship between 
these M responses and the true responses. 


i > IU„ 1 > a oj„2> IU„2> ■ 





This functional relationship most likely varies from system 
to system and also is based upon system operating conditions. 
To deduce what this relationship is for any system would be 
time consuming. Instead, a generalized regression neural net- 
work is implemented to predict the true response. 

A Generalized Regression Neural Network (GRNN) is a 
form of a radial basis neural network suitable for problems 
involving regression, function estimation, or prediction of a 
continuous value, Wasserman, 1993, Alpaydin, 2004. This is 
in contrast to a very similar network, the Probabilistic Neural 
Network, which is used for classification applications. Unlike 
multi-layer perceptrons, which require training, a GRNN is 
constructed from a training set of example inputs and the 
corresponding outputs. The spread of the radial basis function 
is used to compute the bias. The spread, in effect, controls the 
spread of the radial basis function. The example inputs are in 
the form of an mxn matrix where each of the m rows repre- 
sents an observation and each of the n columns represents a 
feature, or parameter. The corresponding example output is 
an mxl vector. The input is said to be n-dimensional. The 
network is divided into 2 layers. The first, or input layer, 
consists of the example inputs. The output layer consists of 
the example outputs. The network generates its output in the 
following manner. The geometric distance is calculated 
between the newly presented input and each of the example 
inputs in the input layer: 
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This produces a vector of length m. Each element of the 
vector is then multiplied by the bias, which is 0.8326/Spread. 
The vector then goes through the radial basis function. 


The radial basis function produces an output that gets closer 
to 1 as the input gets closer to 0. The resulting vector is a 
ranking of how close each of the example cases are to the new 
input. The normalized dot product is then calculated between 
the vector and the example output vector. This is the network 1 5 
output. 


> = -Yx m 


The NNESD is based upon the GRNN approach. Analyti- 25 
cal testing will provide a preliminary validation of the con- 
cept. 

The NNESD concept was analytically tested on a 2 nd order 
Butterworth filter using the MATLAB® matrix calculating 
computer software BUTTER function and the INL LPM 30 
FreedomCAR Battery Test Manual, 2003, Fenton et al., 2005, 
for the lithium ion battery. An SOS input signal was used for 
both cases. For the filter, the component frequencies were 
varied for each run in order to build up a training set with 10 
component lines to be detected. The filter was run 100 times 35 
and the output and taiget response was calculated. The data 
was randomized and half of the data was used to construct the 
GRNN and the other half was used to verify the network. The 
mean squared error (MSE) of the predicted value was calcu- 
lated and the process was repeated 20 times, shuffling the data 40 
each time. The results are shown in FIG. 21 . The MSE of both 
the standard synchronous detection and the CSD are also 
shown for comparison. The MATLAB® matrix calculation 
computer software code used for this analysis is disclosed by 
Morrison, 2005 , Intelligent Self-Evolving Prognostic Fusion. 45 

The technique was then tested on the LPM FreedomCAR 
Battery Test Manual, 2003, Fenton et al., 2005. The input 
parameters to the model were nominally set as per Table 2 and 
randomly varied by up to 5% each run for 1 00 runs to generate 
the dataset. The number of lines, frequency spread and time 50 
step were held as per Table 2. The same training and testing 
scheme that was outlined above was used. The results are 
shown in FIG. 22. The MATLAB® matrix calculation com- 
puter software code used for this analysis is given in Morri- 
son, 2005, Intelligent Self-Evolving Prognostic Fusion. 55 

This limited analytical validation has shown that the 
NNESD concept will significantly reduce the error in the 
estimate of the frequency components of a given system 
response signal. This concept allows a parallel implementa- 
tion for swept frequency measurements to be made by utiliz- 60 
ing a composite signal of a single time record that greatly 
reduces testing time without a significant loss of accuracy. 
Conclusions 

The physical validation of the CSD or NNESD concept 
will rely heavily on work performed by W. Albrecht in his 65 
thesis research: “Battery Complex Impedance Identification 
with Random Signal Techniques,” Albrecht, 2005. In this 


approach, a National Instruments data acquisition and pro- 
cessing system was used along with a custom analog condi- 
tioning system. The CSD or NNESD algorithm could be 
installed directly on the custom analog conditioning system. 
The system software will be CSD or NNESD rather than 
Noise Identification of Battery Impedance (NIBI). The NIBI 
approach would acquire about 100 time records of the battery 
response to noise current . Clearly, a time record would have to 
be of a length of multiple periods of the lowest frequency of 
interest. The CSD or NNESD approach will acquire one time 
record of a length of multiple periods (exact number is still to 
be determined) of the lowest frequency of interest. The exci- 
tation signal would be generated analytically with software, 
as in the NIBI system. The analytical signal would be precon- 
ditioned with a digital low-pass filter as in the NIBI system. 
The CSD or NNESD system may require an analog filter prior 
to the current driver and after the D/A. This analog filter at the 
current driver could serve as the prime anti-aliasing filter. 
This system will use the same bias compensation approach to 
remove most of the DC battery voltage from the acquired 
signal. Improved noise rejection and increased sensitivity 
could be achieved if the voltage sensing were upgraded to full 
differential via a 4-wire system rather than the 2-wire single- 
ended system of the NIBI. An increase in sensitivity will 
enable a reduction in the level of the excitation signal 
required. It is anticipated that the sampled voltage will be 
processed directly with the CSD or NNESD algorithm. It is 
also anticipated that a system calibration would be done 
exactly as the NIBI system by measuring the impedance of 
the test leads to determine any system measurement offset 
and phase shift. 

Currently, calibrating the CSD system for magnitude 
response is done using current shunts in place of a test battery. 
Multiple shunts encompassing range of response are mea- 
sured and the resulting data are used to generate a linear 
regression calibration curve. This technique is also essen- 
tially a one-point phase calibration with the shunt, by defini- 
tion, having zero degrees of phase shift over the frequency 
range of interest. Since the entire measurement system is a 
series process, phase shift in the excitation signal, the test 
object, the detection amplifier or the processing algorithm all 
look like phase shift. Thus, during shunt calibration, addi- 
tional multi-point phase calibrations can be obtained by intro- 
ducing a calibrated phase shift into each frequency of the SOS 
excitation signal and the response processed in the CSD algo- 
rithm by using the non-phase shifted SOS signal. A multi- 
point linear regression magnitude calibration, as well as a 
multi-point linear regression phase calibration, is obtained. 

The principal attribute of the CSD method is the ability to 
obtain a limited resolution frequency spectrum of a system 
function in real time. Resolution of just a dozen frequency 
lines is traded off to obtain the very desirable feature of 
real-time response. The CSD system is able to reduce cross- 
talk between adjacent frequency lines over and above what 
would occur with only synchronous detection. Nevertheless, 
the number of frequency lines in the single time record must 
be limited. In a particularly preferred embodiment, a user can 
choose to increase resolution at the expense of increased 
response time with either operation. The limited number of 
frequency lines in the SOS is adhered to in a single time 
record. However, for either operation additional time records 
are run with each frequency shifted by a shift factor and the 
additional spectra acquired are interleaved to obtain an over- 
all higher resolution spectrum. Thus, a user can make a cus- 
tom trade-off of spectrum resolution against response time for 
an application. 
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The CSD system processes the system response of the SOS 
excitation via repeated synchronous detection acting on a 
residual signal whereby cross-talk contributors have been 
subtracted out. An alternative embodiment of the method of 
detecting system function of the subject invention, as shown 5 
in the flow chart of FIG. 30, uses a fast summation algorithm 
to process the system response of the SOS excitation at step 
32. In the Fast Summation Transformation (FST) embodi- 
ment, all the frequencies of the SOS are harmonics by powers 
of two at step 26. Additionally, the sample period is also a 
power of two with all the SOS frequencies at step 32. Instead 10 
of multiplying the acquired time record by the sine and the 
cosine of each frequency, the SOS is simply rectified relative 
to the square wave and the 90-degree shifted square wave of 
the desired frequency. When the samples of that processed 
time record are summed, all the octave-related harmonics, 15 
other than the frequency of interest will always sum to zero. 
The resulting “In Phase” and “Quadrature” sums can be eas- 
ily processed to yield the magnitude and phase shift of the 
desired frequency component at step 34. 

The FST algorithm can be used to estimate a battery imped- 
ance spectrum at step 36. The excitation Sum of Sines (SOS) 
signal formation at steps 26, 28, 30, the methodology to apply 
it as a current signal to a test battery and the capture of the 
battery voltage response time record at step 32 is almost 
identical to the CSD approach. However, the difference is, if 
the SOS excitation time record contains only sine waves 25 
whose frequencies are all related by octave and harmonics, 
and the sample time is also octave and harmonically related, 
then if that captured time record is “rectified” relative to one 
of the SOS frequencies at step 34, when that transformed time 
record is summed, it will contain only battery response infor- 3Q 
mation relative to that frequency. To identify a specific battery 
frequency response in the SOS signal, the response time 
record is square wave rectified with a phase relationship rela- 
tive to a sine wave of that frequency at step 36. Then, all the 
points in the time record of that transformed signal are simply 
totaled up and normalized to the number of periods of that 35 
frequency present in the SOS signal. This result becomes a 
numerical parameter m l . The process is performed again 
except the rectification square wave phase is relative to the cos 
of the frequency of interest. This other result becomes another 
numerical parameter m 2 . For all the other frequencies, except 40 
the one of interest, all the samples of those sine waves of the 
transformed record will always total to zero. The amplitude 
and phase response of that frequency are obtained as per the 
following relationships. Equation 1 represents the sampled 
signal component at a specific frequency that is to be 
detected. The amplitude, V p and phase, (|) are the desired 45 
information. 


K in \ (In \ (In \ 

—n + (p I = VpsinJ — n Icos0 + Vpcosl — n Isin0 


( 1 ) 

50 


In Equation 3, the signal has been rectified relative to a sine 
wave of that frequency and all the sample values are totaled 
up. 


m i =V P 


Z s K v " + *) - Vp Z sin ( v" + *) 


(3) 


In Equation 4, the signal has been rectified relative to the 
cosine wave of that frequency and again the samples are 
totaled up. Observe that rectification simply involves chang- 
ing the sign of the sample values relative to the sine wave or 
cosine wave timing. 


(4) 


+ 0 ~Vp ^ + + Vp 2 sil {^ w + ^) 


(5) 


m i = Vp^j ^sin^— njcos0 + cos^— wjsin^j - 

^sin^ — «Jcos0 4- cos^— rtjsin^J 


mi = Vpcos0 




VpsinJ 


zHI^-zHI'j) 


*2 


(6) 


m2 = ^sin^— rejcos^ + cos^^reJsin$>j - 

n = 0 v 

3 W_J 

Vp ^ ^sin^— rejcos^ + cos^— «jsin$>j 3 


Vp ^ ^sin^— «Jcos$> 4- cos^— 


Where: 

V p is the amplitude response of the frequency of interest 
N is the number of samples over a period of the frequency 55 
of interest 

(|) is the phase response of the frequency of interest 
n is the discrete time index 

In Equation 1 , N must be constrained as log 2 (N) must be an 
integer greater than 1 . Additionally, the frequency of interest 6Q 
is given as: 


m 2 = Vpcos0 


zW?"))-|'H?”)) + z,H^)) 


VpsinJ 




J NAr 

Where: At is the sample period. 


Note that the parameters K 1? K 2 , K 3 , K 4 are known for each 
65 frequency and m l5 m 2 are the numerical result of the rectify- 
ing algorithm for each frequency. Then, the magnitude and 
phase at each frequency can be obtained as follows: 
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mi = VpcoscpKi + VpsincpKi [mi ] = K \ , K 2 1[ Vpcos(p ' 

=> 

m 2 = Vpcos<pK 3 + Vpsin^A^ [m 2 J = K 3 , K 4 J[ Vpsin0 


The sine wavefonn of the rectification function with zeros 
(7) is given as: 


V P cos(p 1 = K1K4- K 2 K 3 ’ ^2^3 -K 1 K 4 r mi j 

Vpsintp \ = A/3 ATi [m 2 J 

. K 2 K 3 -KiK 4 ' K 1 K 4 -K 2 K 3 _ 


Vpsi# = Ci Vp^Ci+C 2 

Let: , then: , r> \ 

Vpcosfi = Ci ^ = tan -i[£lj 


N 

1, 0 <n< — 

2 

0, = 

N 

-1, — <n<N 

2 J 


For all other signal frequencies present in the rectification 
process, their values must total to zero whenever the rectifi- 
cation is not their specific frequency. The FST system must 
preserve the property of orthogonality between the different 
frequencies. Additionally, any noise present, even after recti- 
fication, will be mitigated by the summation process. Thus, 2 q 
for this technique, mathematically, cross-talk from adjacent 
frequencies will be zero and the SOS time record length can 
be as short as one period of the lowest frequency. However, 
the FST method only works if the octave harmonic relation- 
ship holds for all frequencies in the SOS including the sample 25 
frequency. This ensures that over a period of any frequency 
present in the SOS there will always be an even number of 
samples. The implementation of the rectification functions 
with discrete signals in a manner that preserves orthogonality 
between different frequencies will be discussed in the next 30 
section. 

The algorithm shown above can be realized as computer 
code. A rectification function is simply a unity amplitude 
square wave. FST uses two forms at each frequency, one with 
the phase relationship of a sine wave and the other with the 35 
phase relationship of a cos wave. In a continuous time 
domain, a perfect unit rectification function makes an instant 
transition from -1 to +1 . In a discrete rectification function, 
there are two means for transitioning from +1 to -1 and the 
reverse: pass through zero at a discrete time step or pass 40 
through zero midway between discrete time steps. In the 
discrete time domain, the implementation of a rectification 
function must preserve FST orthogonality between different 
frequencies (i.e., no cross-talk). 

45 

The sine waveform of the non-zero rectification function is 
given as: 


Where: N is the number of samples over the period. 

The cos waveform of the rectification function with zeros is 
given as: 


1, 0 < n < — 

4 


„ , s N 3 N 

Rc(n) = -1, — <n< — 

4 4 

3 N 

0, n- — — 

4 

3 N 

1, —<n<N 

L 4 J 


Where: N is the number of samples over the period. 

There are two means to implement the non-zero rectifica- 
tion function (Equations 9 and 10), but only one means to 
implement the rectification with zeros (Equations 11 and 12). 
All means give nearly identical results. For the non- zero 
rectification function, the more complicated means that will 
work is to take the average between consecutive sample pairs 
to obtain the samples to be summed. Specifically, if one 
averages between points of the sinusoid going around a circle, 
with the last point of the period averaged with the first point. 
The sign of the non-zero rectification function is applied to 
these averaged samples. As an example, consider an eight 
point unity amplitude discrete sine and cos signal. 

One period of an 8-point discrete sine wave: 


= [0, 0.707, 1.0, 0.707, 0, -0.707, -1.0, -0.7 


One period of an 8 -point discrete cos wave: 


1, 0<«<- 


F 1 - 2 * n<N \ 


= [1.0, 0.707, 0, -0.707, -1.0, -0.707, 0, 0.7 


Where: N is the number of samples over the period. 

The cos waveform of the rectification function is given as: 


Rc(n)= -1, 


N 

0 < n < — 
4 


3 N 

— <n<N 
4 


Where: N is the number of samples over the period. 


Thus, with the sine, the 2-point averages rectified by Equa- 
55 tion 9, the non-zero sine rectification: 

:FST^=[(0+0.707)/2+(0.707+l)/2+(l+0.707)/2+ 

(0.707+0)/2-(0-0.707)/2-(-0.707-l)/2-(-l- 

0.707)/2-(-0.707+0)/2]=2(0.707)+2(1.707) 

Now, the cos rectified by Equation 9, the sine: 

60 FSTcs=[(l+0.707)/2+(0.707+0)/2+(0-0.707)/2+(- 

0.707- l)/2-(-l-0.707)/2-(-0.707+0)/2-(0+ 
0.707)/2-(0.707+l)/2]=0 

Now, with the sine 2-point averages rectified by Equation 
10, the non-zero cos rectification: 

:FSTs , c=[(0+0.707)/2+(0.707+l)/2-(l+0.707)/2- 

(0.707+0)/2-(0-0.707)/2-(-0.707-l)/2+(-l- 

0.707)/2+(-0.707+0)/2]=0 
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Finally, the cos rectified by Equation 10, the cos: 

FSTcs=[(l+0.707)/2+(0.707+0)/2-(0-0.707)/2-(- 
0.707- 1)/2-(- 1-0.707)/2-(-0.707+0)/2+(0+ 
0.707)/2+(0.707+l)/2]=2(0.707)+2(1.707) 

Observe the existence of orthogonality with the sine recti- 5 
fied by the cos, or the cos rectified by the sine. This is an 
interesting but unnecessary result that is not needed by the 
FST algorithm. 

The second means of rectification for the non-zero rectifi- 
cation function, and the preferred, is the simplest for the FST 
technique. For rectification, one just changes the signs of the 
captured time record as per Equations 9 and 10. Consider the 
same eight sample discrete sine and cos waves rectified by the 
second means. First, the sine rectified by the sine, Equation 9: 
rectification: 

:FST^=[0+0.707+l+0.707+0-(-0.707)-(-l)-(- 

0.707)]=4(0.707)+2 

Next, the cos rectified by the sine, Equation 9: 

20 

FSTcs=[l+0.707+0-0.707-(-l)-(-0.707)-(0)- 

(0.707)]=2 

Now, the sine rectified by the cos, Equation 10 

FSTsc=[0+0.707-(l)-0.707-(0)-(-0.707)+(-l)+(- 

0.707)]=- 2 25 

Finally, the cos rectified by the cos, Equation 10: 

FSTcc=[l+0.707-(0)-(-0.707)-(-l)-(-0.707)+0+ 

0.707]=4(0.707)+2 

For the rectification with zeros, Equations 1 1 and 12, con- 30 
sider the same eight sample discrete sine and cos waves 
rectified by a second mean. First, the sine rectified by the sine, 
Equation 1 1 : rectification: 

:FST^=[0+0.707+l+(0.707)+0-(-0.707)-(-l)-(- 

0.707)]=4(0.707)+2 35 

Next, the cos rectified by the sine, Equation 11: 

FST5C=[0+0.707+0-(0.707)+0-(-0.707)+0-0.707]=0 

Now, the sine rectified by the cos, Equation 12 40 

FST5C=[0+0.707+0-0.707-0-(-0.707)-(0)- 

(0.707)]=0 

Finally, the cos rectified by the cos, Equation 12: 

FSTcc=[l+0.707+0-(-0.707)-(-l)-(-0.707)+0+ 45 

0.707]=4(0.707)+2 

Again, observe the existence of orthogonality with the sine 
rectified by the cos or the cos rectified by the sine. This is an 
interesting but unnecessary result that is not needed by the 
FST algorithm. 50 

Example 5 

Analytical Testing of Non-Zero Rectification 

(Preferred Non- Averaging) 55 

The non-zero sine and cos rectification functions in the 
preferred non-averaging means, Equations 9 and 10, respec- 
tively, when implemented into the FST Equations 1 through 
8, preserve orthogonality between frequencies. This is shown 60 
via a MATLAB® matrix calculation computer software code 
that builds a discrete time record of an SOS where the fre- 
quencies are octave harmonics including the sample fre- 
quency. All the sinusoids in the SOS start at time zero, have 
zero phase shift and unity amplitude. The process starts with 65 
all waves being sine waves and being rectified with the sine 
rectification function, Equation 9, for each frequency and 
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then repeat with rectification by the cos rectification function, 
Equation 10. Then, the SOS is all cos waves and is rectified 
with the sine rectification function, Equation 9. Cross-talk is 
always expected to be zero and, thus, there is orthogonality 
between frequencies. For this evaluation, the SOS is picked to 
consist of 10 distinct frequencies, all of which are octave and 
harmonically related with the highest frequency having 32 
samples per period. The time record length is 1 period of the 
lowest frequency. For Case 1, the SOS is made from discrete 
sine waves and for Case 2, the SOS is made from discrete cos 
waves. Orthogonality is proven by deleting an arbitrary fre- 
quency from the SOS or SOC and then trying to detect it with 
the FST Equations 1 through 8. If the result of detection is 
zero, then there is no cross-talk from any of the other frequen- 
cies and orthogonality is established. 

Case 1 

SOS is unity amplitude sine waves, 10 octave frequencies 
with 32 samples per period of the highest frequency. The SOS 
time record length is 1 period of the lowest frequency. 
Orthogonality is tested between frequencies by deleting from 
the SOS an arbitrary frequency and then using FST to detect 
it. Anything that is detected will be cross-talk corruption. The 
3rd frequency was deleted. The results for sine and cos rec- 
tification FST are: 

Outs=[0.6366 0.6366 0.0000 0.6366 0.6366 0.6366 0.6366 
0.6365 0.6361 0.6346] 

0utc=[-0.0001 -0.0002 0.0000 -0.0010 -0.0020 -0.0039 
-0.0078 -0.0156 -0.0313 -0.0625] 

As expected, there is no response at the 3rd frequency. In 
fact, whichever frequency is deleted, there will be no response 
at that frequency; thus, orthogonality between frequencies is 
considered valid relative to an SOS. 

Case 2 

The parameters are the same as in Case 1, except SOS is 
cosine waves (SOC). Now, if an arbitrary frequency is deleted 
from the SOC, detection of the deleted frequency measures 
cross-talk. If there is orthogonality between frequencies, the 
result should be zero. The 5th frequency was deleted and the 
cos and sine results, respectively, are: 

Outc=[0.6366 0.6366 0.6366 0.6366 0.0000 0.6366 0.6366 
0.6365 0.6361 0.6346] 

0uts=[0.0001 0.0002 0.0005 0.0010 0.0000 0.0039 0.0078 
0.0156 0.0313 0.0625] 

Again, as expected, there is no response at the 5th fre- 
quency. In fact, whichever frequency is deleted, there will be 
no response at that frequency; thus, orthogonality between 
frequencies is considered valid relative to a SOC. 

Example 6 

Non-Zero Rectification FST Applied to a Battery 
Model 

The method of the subject invention using the FST algo- 
rithm with non-zero rectification was applied to a battery 
model. Specifically, by recursive implementation of the INL 
Lumped Parameter Model (LPM) (FreedomCAR Battery 
Test Manual, 2003). LPM parameters were the same as those 
used in Example 3. Several cases were analyzed and results 
plotted as FST compared to ideal LPM impedance (classical 
joo circuit analysis). The plots were Bode (magnitude and 
phase) and Nyquist. All these plots are shown in continuous 
format rather than the discrete format of Example 3. Never- 
theless, the individual points as computed are clearly indi- 
cated. 
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Case 1 

In the first case analyzed, the SOS started at 0.01 Hz, had 1 0 
frequencies, the time record length was 1 period of the 0.01 
Hz, and the highest frequency was set to 4 samples as this is 
considered a worst-case lower limit. FIGS. 24A-24C show 5 
the impedance spectrum results, respectively, as magnitude 
vs. frequency, phase vs. frequency and the Nyquist plot of 
imaginary component vs. real component. 

Case 2 

It is suspected that the observed error between FST and the to 
ideal response is due to the transient response of the LPM. 
Thus, this case is the same as Case 1 except the time record 
was increased to 3 periods of the 0.01 Hz frequency. FIGS. 
25A-25C show the results in the same format as Case 1 . 

Case 3 15 

In Case 2, it was observed that with the time record of Case 
1 expanded to 3 periods of the lowest frequency, the error was 
greatly reduced but still present. Case 3 is the same as Case 1 
except the time record is set at 2 periods of the lowest fre- 
quency and then the first period is ignored and only the second 20 
period is processed. The objective was to delete from the time 
record the LPM-corrupting transient response. As seen in the 
following plots, FIGS. 26A-26C in the same format as Case 1 , 
the approach worked as the FST and ideal results overlay each 
other. This observation is very strong evidence that the LPM 25 
transient response is corrupting the non-zero rectification 
FST algorithm. 

Example 7 

30 

Analytical Testing of the Rectification Function with 
Zeros 

The second option of a discrete rectification function for 
transitioning from +1 to -1, and the reverse, and passing 35 
through zero at a discrete time step is examined. Equations 1 1 
and 12, when implemented into the FST algorithm and Equa- 
tions 1-8, are shown to preserve orthogonality between fre- 
quencies. The sine and cosine waveform of such a function 
are given, respectively, as Equations 1 1 and 12. These sine 40 
and cosine rectification functions with zeros preserve 
orthogonality between frequencies shown via a MATLAB® 
matrix calculator computer software code that builds a dis- 
crete time record of an SOS where the frequencies are octave 
harmonics including the sample frequency. All the sinusoids 45 
in the SOS start at time zero, have zero phase shift and unity 
amplitude. All sine waves are used initially and are rectified 
with the sine and cosine rectification functions for each fre- 
quency, as shown in Equations 1 1 and 12. Then, the SOS is all 
cosine waves, SOC and rectified with the sine and cosine 50 
rectification functions, as shown in Equations 1 1 and 12. The 
cross-talk results are expected to always be zero. For this 
evaluation, the SOS is picked to consist of 10 distinct frequen- 
cies, all of which are octave and harmonically related with the 
highest frequency having 32 samples per period. The time 55 
record length is 1 period of the lowest frequency. For Case 1, 
the SOS is made from discrete sine waves and, for Case 2, the 
SOS is made from discrete cos waves. Orthogonality is 
proven by deleting an arbitrary frequency from the SOS or 
SOC and then trying to detect it with the FST Equations 1 60 
through 8. If the result of detection is zero, then there is no 
cross-talk from any of the other frequencies and orthogonal- 
ity is established. 

Case 1 

Now, if an arbitrary frequency is deleted from the SOS, 65 
cross-talk is measured when trying to detect the arbitrary 
frequency. If there is orthogonality between frequencies, the 
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result should be zero. The 2nd frequency was deleted and the 
sine and cos results respectively are: 

Outs=[0.6366 0.0000 0.6366 0.6366 0.6366 0.6366 0.6366 
0.6365 0.6361 0.6346] 

Outc=l .0e-01 3 *[0.0003 0.0028 0.0003 -0.0074 -0.0132 
-0.0243 -0.0501 -0.1000 -0.1994 -0.3983] 

As expected, there was no response at the 2nd frequency or 
any other frequency that might be deleted from the SOS, thus, 
there is orthogonality between frequencies. Additionally, 
there are 13 orders of magnitude difference between sine and 
cos rectification results; thus, orthogonality also exists 
between sine and cosine for an SOS. This feature is of interest 
but not necessary for the subject method employing the FST 
algorithm. 

Case 2 

Parameters are the same as in Case 1 , except SOS is cosine 
waves (SOC). Deleting an arbitrary frequency from the SOC 
and then trying to detect it measures cross-talk. If there is 
orthogonality between frequencies, the result should be zero. 
The 7th frequency was deleted and the cos and sine results, 
respectively, are: 

Outc=[0.6366 0.6366 0.6366 0.6366 0.6366 0.6366 -0.0000 
0.6365 0.6361 0.6346] 

Outs=1.0e-01 3 *[-0.0049 -0.0277 0.0051 0.0140 0.0129 
0.0234 0.0032 0.1002 0.1992 0.3983] 

As expected, there is no response at the 7th frequency or 
any other frequency that might be deleted from the SOC. 
Thus, orthogonality also holds for the SOC. Additionally, 
there are 13 orders of magnitude difference between cos and 
sine rectification results; thus, orthogonality between sine 
and cos for an SOC also exists, an interesting but unnecessary 
feature. 

Example 8 

Rectification with Zeros FST Applied to a Battery 
Model 

The rectification with zeros version of the FST embodi- 
ment was evaluated with the recursive implementation of the 
INL LPM exactly as in Example 6. 

Case 1 

In the first case analyzed, the SOS started at 0.01 Hz, had 10 
frequencies, the time record length was 1 period of the 0.01 
Hz, and the highest frequency had 1 6 samples. Results shown 
in FIGS. 27A-27C show the impedance spectrum results, 
respectively, as magnitude vs. frequency, phase vs. frequency 
and the Nyquist plot of imaginary component vs. real com- 
ponent. 

Case 2 

Case 2 was likewise run with the parameters of Example 6 
(the time record is 3 periods of the lowest frequency). Results 
are shown in FIGS. 28A-28C in the same format as Case 1 . 
Case 3 

Again, Case 2 for rectification with zeros is virtually the 
same as non-zero rectification. Corruption by LPM transient 
response is suspected. Case 3 was run exactly as Case 3 for 
non-zero rectification (a time record of 2 periods of lowest 
frequency and process with only the second period via zero- 
rectification FST). The results with the FST algorithm that 
uses rectification functions with zeros is virtually identical to 
the non-zero rectification-based FST in that eliminating the 
first half of the 2-period time record totally eliminates the 
corruption of results by the LPM transient response. Results 
are shown in FIGS. 29A-29C in the same format as Case 1 . 
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Conclusions 

Zero and non-zero rectification for the subject method 
using the FST system have both been validated analytically. 
The rectification without zeros using the preferred means of 
non-averaging is probably easier to implement for FST as it 5 
involves preprocessing the time record by just sign changes, 
and then just summing all the samples. 

It is to be understood that the foregoing examples are 
merely illustrative of the present invention. Certain modifi- 
cations of the articles and/or methods employed may be made to 
and still achieve the objectives of the invention. Such modi- 
fications are contemplated as within the scope of the claimed 
invention. 
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to the frequencies within the sum of sinusoids and at 
least four times the highest frequency in the sum of 
sinusoids; 

rectifying the response time record relative to a sine at each 
frequency of the octave harmonic frequencies with a first 
sine rectification sum m ll5 comprising: 

N i , 

~T~ l N l~ l 

mu = sum = £ TR(nAt) - y TR{nAt) 

n = 0 Ni 

n =~2~ 
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The invention claimed is: 

1. A method for detecting function of a unit under test by 

measuring frequency response the method comprising; the 

acts of: 50 

setting measurement conditions comprising choosing a 
Root Mean Square (RMS) level of current, identifying a 
start frequency, and choosing a period of a lowest fre- 
quency, 

selecting a number of octave harmonic frequencies relative 55 
to the start frequency over which the function of the unit 
under test will be tested, wherein the number of octave 
harmonic frequencies is one or more; 
assembling an excitation time record including a sum of 
sinusoids of the number of octave harmonic frequencies 60 
and a duration of the period of the lowest frequency; 
conditioning the excitation time record to be compatible 
with the measurement conditions; 
exciting the unit under test with the excitation time record 
and simultaneously capturing a resulting response time 65 
record with a data acquisition system having a sample 
frequency that is octave related and harmonically related 


Where: TR is a captured response time record with discrete 
time index n and duration of one period of a first fre- 
quency fi 

N. is a number of sample points for a period of the first 
frequency f : wherein for each sum of sinusoid frequen- 
cies f z , N z - is constrained as log 2 (N z ) and is an integer 
greater than 1 and wherein the frequency of interest is 
given as: 

1 


wherein At is a sample period 

summing the sine rectification sum for each frequency of 
interest of the octave harmonic frequencies to create a 
sine-rectified response time record sum, normalizing the 
sine-rectified response time record sum to a number of 
periods known to be within the sine-rectified response 
time record sum for each frequency, and storing a nor- 
malized sine sum result, m l z -, where i is an \ th frequency 
of interest; 

rectifying the response time record relative to a cosine at 
each frequency of the octave harmonic frequencies with 
a first cosine rectification sum m 21 , comprising: 


N 1 , 3A6 . 

■4 1 N i~ l 

m 2 i=sum= ^ 77?(reAr) - ^ TR(nAt) + ^ TR(nAt) 

n = 0 /V, 3 Ni 

n =- 4 n =~r 


summing the cosine rectification sum for each frequency of 
interest of the octave harmonic frequencies to create a 
cosine-rectified response time record sum, normalizing 
the cosine-rectified response time record sum to the 
number of periods known to be within the cosine-recti- 
fied response time record sum for each frequency, and 
storing a normalized cosine sum result, m 2 i , where i is 
the \ th frequency of interest; 

determining a real and an imaginary form of a magnitude 
and a phase of each of the octave harmonic frequencies: 


' V/cos^i ' 

Vising,- 


*4 

KiK a -K 2 K^ 

k 3 

K 2 K 3 -KiK 4 ' 


k 2 

K 2 K 3 -KiK 4 

Ki 

KiK 4 -K 2 K 3 . 


[ m u 
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Where: 

V z is a magnitude of the voltage for the frequency f z 
(|) z is a phase of the voltage for the frequency f z 

m 2 . are rectification sums for the frequency f. 5 


-s'Ka-SWr)) 


Ni 

”=T 


10 


Rc{n) - 


1, 0 < n < — 

4 


N 3 N 

- 1 , — <«< — 
4 4 

3 N 

0, „ = - 

3 N 

1 —— <n < N 

4 




K 2 


= SHI-))- SHI-)) 

n= Y 


15 


¥-1 


3^-1 


^=zKr))-2Kr)) + zWI")) 

m . w. 

n =H 


N i 

n =i 
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,= S Hr))- £ Hr)) + % Hr)) 

,,=0 


Ni~ 1 
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and, obtaining the magnitude and the phase via a standard 
rectangular to polar conversion: 35 


V;sin0; = Cl 

Let: , then: 

V/COS0; = C2 


V/ = V C 1 +C 2 


pi = tan' 


ID 


40 


4. The method of claim 1, further comprising the acts of: 

repeating the original acts of; selecting the number of 

octave harmonic frequencies, assembling the excitation 
time record, conditioning the excitation time record, 
exciting the unit under test, rectifying the response time 
record relative to the sine at each frequency, summing 
the sine rectification sum for each frequency of interest, 
rectifying the response time record relative to the cosine 
at each frequency, summing the cosine rectification sum 
for each frequency of interest, determining the real and 
the imaginary form, and assembling the magnitude and 
the phase, at least once, wherein the octave harmonic 
frequencies are shifted in the repeated acts relative to the 
original acts; and 

interleaving the frequency response of the original acts 
with the frequency responses of the repeated acts to 
obtain a higher resolution frequency response. 

5. The method of claim 1, wherein: 

rectifying the response time record relative to a sine at each 
frequency of the octave harmonic frequencies includes a 
second sine rectification sum m 12 , comprising; 


sum = 


AH , N 2 , 

~ 2 l ~ 1 ^2—1 /V 2 +-f-l 2/V 2 -l 

^ TR(nAt) - ^ TR(nAt) + TR(nAt) - TR(nAt) 

n= 0 N 2 n=N 2 .. N 2 

n= ~T L n=N 2 + ~T 


assembling the magnitude and the phase at each of the 
octave harmonic frequencies to determine the frequency 
response. 

2. The method of claim 1, wherein rectifying the response 
time record relative to the sine at each of the octave harmonic 
frequencies is performed for signals that pass through zero at 
a discrete time step, wherein N is a number of samples over 50 
the period: 


sum 
m 12 =— ; 

and 

rectifying the response time record relative to a cosine at 
each frequency of the octave harmonic frequencies 
includes a second cosine rectification sum m 22 , compris- 
ing; 


Rs(n) = 


0, n = 0 

N 

1, 0 <n< — 

2 



- 1 , 


N 

2 


<n < N 


55 


60 



sum= ^ 

n= 0 


3 /V 2 , 

-4 2 - 1 "2-1 

TR(nAt) - TR(nAt) + TR(nAt) 4- 

N 2 3AH 

*=-r n= ~r 

A^+^f-l N 2 +-j L ~l 2A/ 2 -l 

^ TR{nAt) - Yj TR(nAt) + ^ TR(nAt) 

n=N 2 . 1 AH 3N 2 

2 n =N 2 +-f- n =N 2 +- 4 ± 


3. The method of claim 1, wherein rectifying the response 
time record relative to the cosine at each of the octave har- 
monic frequencies is performed for signals that pass through 65 
zero at a discrete time step, wherein N is a number of samples 
over the period: 


sum 
m 22 = — ; 

where: N 2 is a number of sample points in a second fre- 
quency at an octave harmonic of the first frequency and 
2N 2 =N 1 . 
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6. The method of claim 5, wherein: 
rectifying the response time record relative to a sine at each 
frequency of the octave hannonic frequencies includes a 
third sine rectification sum m 13 , comprising; 



sum= y 

12=0 


/V 3 -l 

TR(nAt) - TR(nAt) + 

N-x 

n =2 


N-x Nx 

N 3 +-^—l 2 /V 3 -I 2N 3 +-j--i 

2 TR(nAt) - ^ TR(nAr) + £ TR(nAt) - 

n=N 3 N 3 n=2N 3 

5 n=N 3 +--f- J 


Nx 

3/V3-I 3N 3 + -f-l 4N 3 -l 

TR(nAt) + ^ TR(nAt) - ^ TR(nAt) 

N 3 n=3N 3 

n=2N 3 +-ig- * n= 3 A ^3 + --y- 


sum 
m i3 = — ; 


and 

rectifying the response time record relative to a cosine at 
each frequency of the octave harmonic frequencies 
includes a third cosine rectification sum m 23 , compris- 
ing, comprising: 


sum = Z z 

12=0 n _N^ 


1 1 - ^ 

3Nx 


3N 3 

n =~T 


setting measurement conditions comprising choosing a 
Root Mean Square (RMS) level of current, identifying a 
start frequency, and choosing a period of a lowest fre- 
quency; 

5 determining a theoretical frequency response for the refer- 
ence standard; 

selecting a number of frequencies and a frequency spread 
over which the reference standard will be tested, 
wherein the frequency spread is by octaves; 

10 assembling a phase-shifted excitation time record includ- 
ing a sum of sinusoids of the number of frequencies with 
a pre-selected phase shift, wherein a duration of the 
excitation time record is greater than or equal to one 
15 period of a lowest of the number of frequencies; 

conditioning the phase-shifted excitation time record to be 
compatible with the measurement conditions; 
exciting the reference standard with the phase-shifted exci- 
tation time record and simultaneously capturing a 
20 response time record with a data acquisition system 
having an appropriate sample frequency; 
processing the response time record to obtain estimated 
magnitude and estimated phase at each of the number of 
frequencies to determine the frequency response; 

25 wherein, phase calibration is achieved through a selection 
of multiple known phase-shifts until enough points are 
obtained to determine phase calibration constants for 
each of the number of frequencies using a regression 
analysis. 

1 1 9. The method of claim 8, wherein the act of processing the 

response time record includes: 

rectifying the response time record relative to a sine at each 
of the number of frequencies, with a first sine rectifica- 
„ tion sum m l l5 comprising: 


/v 3 +^p--l N 3 +--£--\ 2 /V 3 -I 

^ TR(nAt) - Yj TR(nAt) + ^ TR{nAt) + 

n=N 3 N 3 3N 3 

3 n=N 3 +-f n=N 3 + -jf- 

2N 3 ^--i 2 W 3 + - 3 ^ 2 .-i 3^-1 

2 TR(nAt) - ^ TR(nAt) + ^ TR(nAt) + 

n=2N 3 ... N 3 ... 3N 3 

* n=2N 3 +-£- n=2N 3 + 

3N 3 + ^--l 3/V 3 + 3 -^-l 4 /V 3 -I 

Yj TR(nAt) - Yj TR{nAt) + ^ TR{nAt) 

n=3N 3 xu ^ N 3 xu J^3 

n=3N 3 + -f n=3N 3 + ~^ 

sum 

m 23 = — ; 

where: N 3 is a number of sample points in a third frequency 
at an octave harmonic of the first frequency and 4N 3 =N 1 . 

7. The method of claim 6, comprising: 
rectifying the response time record relative to a sine at each 

frequency of the octave harmonic frequencies includes 
at least one additional sine rectification sum m lw ; and 
rectifying the response time record relative to a cosine at 
each frequency of the octave harmonic frequencies 60 
includes at least one additional cosine rectification sum 
m 2 „; wherein a last frequency includes 2”' 1 periods 
within a period of the first frequency. 

8. A method of calibrating phase response of a captured 
time record for a unit under test by measuring frequency 65 
response of a reference standard the method comprising the 
acts of: 


1 Wi-i 

mu = sum = y TR(nAt) - y TR(nAt) 

Art n=0 Nt 

40 n=- f 

Where: TR is a captured voltage time record with discrete 
time index n and duration of one period of a first fre- 
45 quency f : 

Nj is a number of sample points for a first frequency f l 
wherein for each sum of sinusoid frequencies f z , N z is 
constrained as log 2 (N f ) and is an integer greater than 1 
and wherein the frequency of interest is given as: 

50 

* = N^At 

wherein At is a sample period 

summing the sine rectification sum for each frequency of 
interest of the number of frequencies to create a sine- 
rectified response time record sum, normalizing the 
sine-rectified response time record sum to a number of 
periods known to be within the sine-rectified response 
time record sum for each frequency, and storing a nor- 
malized sine sum result, m ljf , where i is the \ th frequency 
of interest; 

rectifying the response time record relative to the cosine at 
each of the number of frequencies, with a first cosine 
rectification sum m 21 , comprising: 
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, 3/V, , 

-1 -4^-1 N { - 1 

r TR{nAt)~ Y TR(nAt) + Y TR(nAt) 

=o /Vi 3/Vi 

"=X n =4 


10. The method of claim 9, wherein: 

rectifying the response time record relative to a sine at each 
of the number of frequencies includes a second sine 
rectification sum m 12 , comprising: 


summing the cosine rectification sum for each of the num- 
ber of frequencies to create a cosine-rectified response 
time record sum, normalizing the cosine-rectified 10 
response time record sum to the number of periods 
known to be within the cosine-rectified time record sum 
for each frequency, and storing a normalized cosine sum 
result, m 2 z , where i is an \ th frequency of interest; 

determining the real and imaginary form of the magnitude 
and phase of each octave harmonic frequency: 


sum = ^ TR(nAt) - TR(nAt) + Y TR(nAt) - Y TR(nAt) 


Vicoscpi 1 k { k 4 - k 2 k 3 ’ k 2 k 3 -k { k 4 [ m,\ 

V/sin^; J K 3 Ki [ m 2 j 

_ K 2 K 3 -K x K 4 ' K[ K 4 — K 2 K 3 . 


rectifying the response time record relative to a cosine at 
each of the number of frequencies includes a second 
cosine rectification sum m 22 , comprising: 


V z is a magnitude of the voltage for the frequency f z 
(|) . is a phase of the voltage for the frequency f. 
m 1 ., m 2 . are rectification sums for the frequency f. 


sum = Y TR(nAt) - Yj TR(nAt) + Y TR(nAt) + 

n=0 _N 2l _ 3 N 2 


Y TR(nAt) - Y TR(nAt) + Y TR{nAt) 


K ' = z HS"))~ z KH) 


^£Hr))-SHr)) 


^=zKI«))-T.K^)) + SHr)) 


where: N 2 is a number of sample points in a second fre- 
40 quency at an octave harmonic of the first frequency and 

2N 2 =N 1 . 

11. The method of claim 10, wherein: 
rectifying the response time record relative to a sine at each 
45 of the number of frequencies includes a third sine recti- 
fication sum m 13 , comprising: 


50 sum = Yj TR(nAt) - Y TR(nAt) + 


k <= z HsH) - z HsH) + z HU) 


N 3 +-*—l 2/V 3 -l 2N 3 +-£--l 

Y TR{nAt) - Y TR(nAt)+ Y TR{nAt) - 

n=/V 3 * 3_ n=2N 3 


and, obtaining the magnitude and the phase via a standard 
rectangular to polar conversion: 


Y TR(nAt) + Y TR{nAt) - Y TR{nAt) 


y-sin *= Cl v i=^ c2 i+ c l 

Let: , then: , n N 

V ; cos^ = C 2 0,'=tan-‘(g). 
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and 


rectifying the response time record relative to a cosine at 
each of the number of frequencies includes a third cosine 
rectification sum m 23 , comprising, comprising: 


32 

-continued 

3W3 + -^--1 3/V3 + — ^ — 1 4/V3-I 

£ TR(nAt) - Yj TR(nAt) + ^ TR(nAt) 

n= 3/V3 N 3 3/Vo 

3 n=3/V3 + -^- n= 3 N 3 -\ — 


N 3 , 3/V3 , 

•^■-1 /V 3 -l 

sum = TR(nAt) - ^ TR(nAt) + ^ TR(nAt) + 

n= 0 N 3 3/V3 

A/3+^2.- 1 /V 3 + --^2--l 2/V3-I 

^ TR(nAt) - ^ TR(nAt) + TR(nAt ) + 

«=/V 3 a, n 3 a, 3A, 3 

3 n=N 3 +~f n =N 3 + -^- 

2/V 3 +^--l 2/V 3 + - 3 ^-l 3/V3-I 

2 77?(«Ar)- ^ 77?(«Af) 4 - ^ 77?(«Ar) 4 

n=2N 3 , /V 3 3/V3 

3 n=2N 3 + n=2N 3 + 


10 


15 


where: N 3 is a number of sample points in a third frequency 
at an octave harmonic of the first frequency and 4N 3 =N 1 . 
12. The method of claim 11, wherein: 
rectifying the response time record relative to a sine at each 
of the number of frequencies includes at least one addi- 
tional sine rectification sum m lw ; and 
rectifying the response time record relative to a cosine at 
each of the number of frequencies includes at least one 
additional sine rectification sum m 2w ; 
wherein a last frequency includes 2”' 1 periods within a 
period of the first frequency. 



